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ABSTRACT 


Interconnectors in modem VLSI structures are getting miniaturized in cross- 
sectional area and in length. During processing of these structures, residual 
stresses are generated in the aluminium interconnections due to the thermal 
loading. These stresses are of very high magnitude (hundreds of MPa) and are 
responsible for failure of these aluminium interconnections. Some of the most 
significant failure mechanisms are stress migration, plastic flow and thermal 
grooving. 

In present work, failure analysis of the above structure is being carried 
out with the help of Boundary Integral Method. The problem is being treated in 
semi coupled manner, i.e. for each time step, while cooling, a thermal analysis is 
performed and then the resulting temperature profile over the domain is used for 
stress field estimation. For transient analysis, time dependent fundamental 
solution scheme is used in the BEM. The Von-Mises yield criterion in plasticity 
is applied to examine whether the relaxation of stresses may be caused by plastic 
deformation. Focus has been made for maintaining the integration accuracy at 
sub-micron domains. The stresses-induced failure in an aluminium conductor is 
investigated in terms of diffusion along the surface and grain boundaries. The 
mechanism of thermal grooving, which initiates and grows at elevated 
temperature, is modeled and failure process is numerically simulated. Effect of 
the stress concentration near the groove tip is considered while analyzing 
thermal grooving. 
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Chapter 1 

INTRODUCTION 


LI Introduction: 


Interconnectors in modem very large integrated circuits (VLSI) are getting 
miniaturized wherein the thickness of the conductors is reaching dimensions of 
less then ljim. It, therefore, become more essential for assuring reliability to 
understand the failure mechanisms at submicron level, which differs greatly 
from those of the components of large dimensions. During fabrication, these 
circuits are subjected to several thermal cycles. Thus, residual stresses generate 
in the aluminium interconnections due to the thermal loading. These stresses are 
of very high magnitude (hundreds of MPa) and are responsible for failure of 
these aluminium interconnections. Some of the most significant failure 
mechanisms are plastic flow, thermal grooving and stress migration. 



Fig 1.1: Final multilevel structure 
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It is well known that stress field and high temperature causes the 
diffusion of the atoms (or vacancies) along a grain boundary. This migration 
results in creep cavity and gives rise to a groove at the grain boundaries. The 
grooves generated prorogate even at the low temperatures due to the residual 
stresses in the interconnectors. This process is known as thermal grooving. 



Fig 1.2a: Voids observed in A1 line which causes “Stress Migration” 

Failure 



Fig 1.2b: Photograph of the aluminium conductor in the VLSI package which 
failed due to stress induced migration 

Stress migration is one of the principle mechanisms that cause failure of 
aluminium wire. This problem was first identified by Klema et al.[l] and Curry 
et al.[2]. This failure is caused by the nucleation and growth of the voids in the 
aluminium line under conditions of high stress at elevated temperatures. This 
process does not take into consideration the effect of current flow in the 
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aluminium line, where the electron wind provides a driving force to the 
aluminium ions. This force drives the ions to move in the conducting line and 
cause open-circuit-failure. 

The mechanical stress which causes migration of atoms is induced 
during the line passivation, wherein the metal lines are covered with a glass film 
at high temperature and subsequently cooled down to room temperature. Since 
the thermal expansion coefficient of aluminium is much higher than that of the 
passivation films (approximately 20 times larger), the cooling process generates 
large thermal stresses. The existence of residual stresses is detrimental to the 
reliability of the aluminium interconnectors as it enhances the failure mechanism. 

The boundary element method, based on the boundary integral equation 
formulation is well known for its accuracy and efficiency in stress analysis. In 
recent years, BEM has also been found to be especially accurate and efficient in 
the analysis of thin films. FEM studies show that very fine mesh is essential in 
order to accurately evaluate the interfacial stresses in thin films and coatings. 
For cases involving singular fields, such as interfacial cracks or free edges, even 
finer mesh are required, which can quickly make the domain based FEM 
approach costly. 

Any analysis methods to assess the residual thermal stress patterns in a 
multilevel VLSI structure encounters two difficulties: first, the structure is 
subjected to several thermal cycles due to the various processing steps involved; 
secondly, the geometry of the structure itself changes from one fabrication step 
to another (as material layers are added or removed). 

In the present work an advanced strategy has been applied which 
overcomes all previously mentioned drawbacks. With only one mesh and only 
one run of the program (using a so called activating-deactivating strategy) this 
strategy shows how to obtain the complete stress distribution during the 
manufacturing of the multilevel structures. The developed BEM approach can 
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provide not only a robust numerical tool for the analysis of multilayer systems, 
but also the basis for further investigations of interfacial cracks, thermal effects 
and contact mechanics for the multi-layer system or various thin films. For 
transient thermal analysis, time dependent fundamental scheme has been used. A 
stress induced migration failure of micron or submicron aluminium line in terms 
of interaction between diffusion along the surface and grain boundary takes 
place at high temperature levels. A numerical modeling strategy for estimation 
of such phenomenon has been developed which can be also be used to time to 
failure prediction. 

1.2 Literature Survey: 

Over the past three decades, due to miniaturization of electronic devices, the 
width of the metal interconnection lines used in integrated circuits has decreased 
substantially. The accompanying increase in the line aspect ratio (defined as the 
ratio of layer height and line width) affects the evolution of the thermal stresses, 
which are generated from thermal expansion mismatch between aluminium line 
and the surrounding materials. Numerical modelling using the finite element 
method has been used to predict the evolution of the stresses and the effect of 
line aspect ratio in metal interconnectors (Shen [3]). 

Igic and Mawbay [4] presented an advanced strategy for modeling the 
thermal stresses induced in aluminium interconnectors during processing of 
multilevel structures. He performed the analysis in two steps wherein he used 
the residual stresses from one processing step to be used as the initial condition 
for the next step. Igic and Mawbay [5] presented the effect of titanium barrier 
layer on the residual stress in passivated multilevel aluminium metallization 
structure. 
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Mullins [6] presented a theory which describes the development of 
surface grooves at the grain boundaries of a heated polycrystal. These grooves 
develop on the surface of a hot polycrystal wherever a stationary grain boundary 
intersects the surface due to the diffusion of the atoms through the free surfaces 
and the grain boundaries. These grooves then provide early nuclei for crack 
initialization and propagation through grain boundaries and generate a creep like 
property in the aluminium interconnectors. Martinez and Nix [7] used the 
Mullins theory to develop finite difference technique for the diffusive growth of 
2-D and axisymmetric cavities having equilibrium shapes and located on grain 
boundaries. The shape evaluation and growth kinematics of individual cavities 
as well as the time required for adjacent cavities to grow was formulated as a 
function of applied stress and the ratio of grain boundary to surface diffusivity. 
Kato et al [8] analyzed the relaxation of the residual stresses in aluminium line 
by plastic deformation and diffusion analytically. He extended the Mullins 
theory for the thermal grooving process in aluminium lines. Kitamura, Ohtani 
and Yamanaka [9] investigated the mechanism of grooving which initiate and 
grow at elevated temperature along the grain boundary and numerically 
simulated failure process of aluminium interconnectors. Igic and Mawbay [10] 
suggested an improved numerical model based on finite difference discretization 
for the diffusion relaxation and groove growth in the aluminium line. This 
model is very useful for time to failure prediction and provides better 
understanding of thermal grooving phenomenon. 

In the last three decades, FEM based on plate and shell theories has been 
a successful tool for the analysis of 3-D thin structure such as layered composite 
structures to study their deformation and stress in the micro-scale. However, 
most plate and shell theories are based on various assumptions about the 
geometry, loading and deformation of the structure, and therefore the accuracy 
and reliability of the FEM for the thin structures in the micro or nanoscales are 
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in doubt. The BEM based on the elasticity theory and boundary integral 
equation formulation (Banerjee [11]) is in general more accurate in stress 
analysis of thin structures. 

Pina and Fernandes [12] presented a boundary element formulation for 
three dimensional time-dependent heat conduction problems in homogeneous 
isotropic bodies using a time marching scheme which does not require the 
discretization of domain. Pasquetti, Curso and Wrobel [13] proposed a time 
dependent fundamental scheme to solve transient heat transfer problem. This 
scheme was found very efficient and easy to apply and did not need any domain 
integral. Luo, Liu and Berger [14] developed a model for the application of the 
boundary integral method in the (2-D) thin structures with the thickness to 
length ration in the micron or nanoscales. He presented an efficient analytical 
method to deal with the nearly singular integrals. In addition, he proposed a new 
nonlinear coordinate transformation for nearly singular integrals to further 
increase the numerical accuracies. Luo, Liu and Berger [15] advanced his 
boundary integral model for 2-D thin structure to the analysis of interfacial 
stresses for multi-coating systems. He also supplied a basis for the analysis of 
interfacial cracks, thermal effects and contact mechanics using boundary integral 
method. 

Chen and Liu [16] developed an advanced boundary integral method for 
analyzing thin layered structures under the steady state thermal loading and also 
compared the thermal BEM with commercial FEM software for the multi- 
layered structures of nanoscale dimensions. 

1.3 Present Work and Layout of the Thesis: 


The objective of the present work is to determine the stress field of an 
aluminium layer in multilayer VLSI structure subjected to cyclic thermal 
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loading during and after processing. First the temperature profiles due to various 
fabrication steps have been solved. During thermal loading structure is subjected 
to natural convection boundary conditions. The resulting temperature profile has 
been used to get the residual stresses in the structure. Thermal grooving process 
which occurs due to the residual stresses in aluminium layer has been analyzed. 
Focus has been made for accurate evaluation of domain integrals for sub-micro 
scale domain. The boundary integral method has been made use of for 
accurately performing the above task. 

The thesis has six chapters and each of them is self contained unit. 

In Chapter 2, the 2-D transient heat conduction formulation is described 
using time dependent fundamental scheme. This chapter also includes a 
discussion on the treatment of various types of initial conditions and an efficient 
convolution-type time-marching scheme. 

Chapter 3 deals with the thermal stress analysis for two dimensional 
problems. Thermal effects over the solution domain are treated to calculate 
stress fields. Effect of the multiple domains and special consideration for 
evaluating integrals at the sub-micron scale are explained. A simplified model 
for the stress relaxation due to plastic deformation is also explained. 

In Chapter 4, the stress induced failure in aluminium interconnectors is 
investigated in terms of diffusion along the surface and along the grain boundary. 
A numerical model dealing with thermal grooving process has also been 
described. 

In Chapter 5, the results have been presented for stress analysis during 
processing and thermal grooving phenomenon after the processing. 

Finally, in Chapter 6, the work is summarized and the relevant 
conclusions and the scope for the future work are presented. 
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Chapter 2 

TRANSIENT TEMPERATURE ANALYSIS AND BEM 


This Chapter describes boundary element formulation for transient heat 
conduction using time dependent fundamental solution. Numerical techniques 
are then employed to solve the integral equations in discrete form through time 
marching procedure. Space and time discretization algorithm are discussed in 
detail for two dimensional problems and expressions resulting from the 
analytical integration of the fundamental solutions with respect to time are 
included. The next section describes how linear boundary conditions are 
implemented. 

The BEM formulation is given in brief and all salient points of the 
method relevant to present problem have been brought out clearly. This 
temperature field will be used to evaluate residual stresses and thermal grooving 
in the next chapters. 

2. 1 Problem Definitions: 

The governing differential equation for transient heat conduction can be written 
as [15] 

V 2 r+6 = i— (2.1) 

k & 

where T is temperature, t is time, k is thermal diffusivity (k = X/pc), X is the 
thermal conductivity, p is density, c is specific heat, and the term b represents 
the internal heat generation. 

The most common boundary conditions in heat transfer problem are of 
the following types: 
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• Dirichlet condition (prescribed temperature): 

T = T on n (2.2) 

• Neumann condition (prescribed flux): 

dT 

q= a on r 2 (2.3) 

on 

• Fourier condition (convection): 

p(T) = -h(T-T a ) + t on Ts (2.4) 

in which h is the heat transfer coefficient, T a is the ambient temperature and <D is 
a known value. The adopted convection implies that q is positive if the surfaces 
flux p = Aq is inwards the region. Moreover, a and p may vary in space and time. 



Fig 2.1: A domain V with boundary F 

/ 

The present problem is being solved in two steps. In first step BPSG 
layer is put over the silicon layer at passivation temperature and cooled down to 
room temperature. While in the Step-2 aluminium layer passivated with PSG 
layer is deposited over the BPSG layer. Thus, for the first step, the initial 
temperature profile is constant passivation temperature, while in the next step 
the temperature profile calculated from first step is being used as the initial 
condition. 
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2.2 BEM Introduction: 


The BEM is a numerical technique based on integral equation formulation of the 
continuum mechanics problems. There are two types of integral equation 
formulation. One contains, as basic unknowns quantities with a clear physical 
meaning. This type is called the direct formulation. In the second type, known as 
indirect formulation, the basic unknown quantities have no physical meaning. 
The physical variables are obtained from these quantities. In present work, 
numerical approach used is based on direct formulation [15]. 

There are some general characteristics of the BEM which lead to 
advantages for the analysis of static and dynamic continuum mechanics 
problems. First of all, as the equations are formulated in terms of variables on 
the boundary, only the boundary has to be discretized as opposed to domain 
techniques, such as Finite Element Method (FEM) and the Finite Difference 
Method (FDM). As a consequence, the resulting number of equations is smaller 
in the BEM. This leads to less data preparation and less requirement of computer 
time and storage. Another consequence is that the mesh generation process is 
simpler than in domain type methods. This advantage is important for the 
situation where the geometry changes throughout the solution process, as in 
crack propagation study, processing of multilayer or multi-coating structures, etc. 
A second characteristic of the BEM is that it produces accurate solutions on the 
boundary and in particular at any selected internal point. This feature makes the 
method very appropriate for problem where high accuracy is required. 

Among the disadvantages of BEM, first is the mathematical complexity 
in BE formulation is more (evaluation of Green function). The solution matrix 
resulting from the BE formulation is unsymmetric and frilly populated with 
nonzero coefficients. This means that the entire solution matrix is needed to be 
saved in the computer memory. Second disadvantage is its difficulty in dealing 
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with non-linear material properties. In such cases the integral equations include 
domain integrals and the advantages of a boundary formulation vanish if there 
are many non-linear material zones. 

2.3 Boundary Integral Equation and Fundamental Solution: 


The governing heat conduction equations (2.1) to (2.4) can be recast as an 
integral equation over space and time [13]. Divergence theorem is used to 
convert the governing differential equation into a boundary integral equation as 
follows: 


c u T M tp + £ [kTq'dTdt= £ [kqT'dTdt+ £ ^kbT'dOdt + jTj'dQ 


(2.5) 


where T M is temperature at point M at time tp, c M is a coefficient depending on 

the geometry of T at point M, k is thermal diffusivity, T temperature profile at 
the boundary f. To is initial temperature profile over the domain, b is body force, 


dT dT* — >• 

q = — , q* - - , and n is the unit outward normal vector. 

dn dn 

The fundamental solution T* describes the temperature field due to a unit 
heat source applied at pint M at time to. Many researches have given different 
fundamental solutions of different accuracies (Chandra, Chan and Lim [18], 
Koizumi, Utamura and Kotani [19]). In the present work the fundamental 
solution given by Wrobel and Brebbia £11} is used, which is given by, 


r*(r,r) = 


1 


exp 




(4n*r)* /2 
Id 

(4kT) (s+2) n sl2 


exp 


( 4 k2 


( 2 . 6 ) 


(2.7) 
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with d = -n.r , r being the position vector connecting source and field points, 

T = t F — t and s is the number of spatial dimension of the problem. 

This formulation can also be applied to non-linear problems in which the 
thermal conductivity X is a function of temperature. In this case, it is convenient 
to employ the Kirchhoff Transformation, 

T 

yr=\MT')dT' (2.8) 

T t 

where T r is the reference temperature. For steady state problems, the above 
transformation is sufficient to linearize the governing equation. For transient 
problems, the integration by parts of the equation (2.1) weighted by T* gives rise 
to additional domain integrals involving the derivatives of k. Assuming, that the 
space and time derivatives of k are small quantities and can be neglected, then, 

— = 0 (2.9) 

8t 

An integral equation similar to (2.2) is obtained for the non-linear 
problem [17] 

C MWu,t F cw qdTdt = £ [kpT'dTdt+ [kgT'dQdt + tyj'dQ. 

( 2 . 10 ) 

2. 4 Numerical Solution of the Boundary Integral Equation: 

The numerical solution of the boundary integral equation (2.5) and (2.10) 
requires space and time discretization. Two different time-marching schemes 
can be used in the solution: 

1 . For each time step, to calculate temperature at the resolution time t F , the 
temperature at the previous the time level t F -i is taken as initial 



temperature. This approach minimizes the time integrations but requires 
the evaluation of the domain integral associated with the temperature 
field Tf-i at each time step; 

2. For each time step, the solution is restarted from the initial time to; in this 
way, the domain integral associated to the initial conditions can be 
avoided for the majority of practical situations. 

The second approach was adopted in this work. Recently, truncation 
algorithms were developed [13] to improve the computer efficiency of such an 
approach. One such algorithm, named ‘steady temporal domain’, is herein 
presented in detail. 

2.4.1 Time and space discretization: 


Equation (2.5) represents two domain integrals due to initial conditions and 
internal source. For simplicity, if there is no internal heat generation and that the 
initial temperature is constant, then it is possible to rewrite the problem into an 
equivalent one with zero initial value for temperature difference (T — To). For 
the non-linear case, it is sufficient to take To as reference temperature T r in the 
Kirchhoff transformation. 

With the above simplification, Eq. (2.5) can be written for a point M in 


the portion ofT as 

\r u , f + ( f j [kTqdTdt = £ [kqT'drdt (2.1 1) 

Dividing the boundary T into N boundary elements and the time span 
(t F - to) into F time steps, the following discretized equation is obtained 

Sl i, [I [MT'dvdt ( 2.i2) 

where T* F being the temperature at node i at time t F . 
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Assuming that the variables are linear in each time step, the temperature 
variation on element j between time interval f-\ and/is given by 

T jM±<tt7»+W-'fJ CU3) 

A i f 

with At f= if- tf.j. A similar expression can be written for q. Writing 




(2.14) 



(2.15) 

G K,.f,, ■ 

= 77- f lk(t,-t)T'Srdt 

A tf_ { *•/-! *j 

(2.16) 


[1 

- 

*r-* 

1 

# 

£*■ 

(2.17) 

Muss 

= h, J s,s + jt,A, / 

(2.18) 


in which 8 is the Kronecker delta, the relations Eq. (2.12) can be written as: 


F N 


F N 


2. J ^) = 2E«a w „^+G2 IJjr ^) (2.19) 

/= /=i 7=i 

Writing the above equation at all boundary nodes / (1</<N) using a 
collocation technique, the following system of equations is obtained 

2(//l F , / r / _, + H2, v r / )=2(Gl F , / e / - 1 +G2 F-/ e / ) (2.20) 

/=1 /= 1 

The temperature and flux at each node are known for all time levels 
previous to t F , so that the above equation can be rewritten in the form 

H2 Ff T f =G2 FF Q F +S F 

with 
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s F =t,(m 

FJ T f - 1 +H2 FJ T f )-H\„T F _, + F f i {G\ 

FjQf-l + G2 f f Qf)+Gl F f Qf - i 


/= i 


/=! 


(2.22) 

The details of evaluation of HI, H2, G1 and G2 can be found in 
Appendix ‘A’. 


2.4.2 Shave functions for the element: 

The boundary of the solution domain is divided into a number of connected 
elements. The variation of geometry and variables over element may be constant, 
linear, quadratic, cubic or higher order and the variation may be different for 
both geometry and the variables. Increasing the order of variation of an element 
produces accurate solutions, but the penalty of higher CPU must be paid. These 
variations of geometry and variables over an element are prescribed by defining 
the “Shape functions”. The geometry of a 3-noded quadratic element can be 
described by the coordinates of its three nodes and the shape function as follows: 

= = n ,(i})x 1 +N 2 (ti)x 2 +N i (r/)x 3 

(2.23) 

y{v) = 2X = N i +n 3 {v)y 3 

C=1 

where the shape function can be defined as 

N 2 =(U-j,)(l-r,) (2.24) 

AT 3 =|(l+7) 
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Fig 2.2: An isoperametric quadratic element 


Isoperaimetric elements are being used in the present analysis. Therefore, 
the shape functions will be same to define the solution variables, as follows: 

T(i)- = ix {n)T c = JV, {v)T, +n 2 (v)r 2 +n 3 (v)t 3 


C=l 




( 8T 


a,l +JVlW 


dn J 3 


(2.25) 


2.5 Assembly of the Different Matrices for Multi-layer System: 

Fig (2.3) and (2.4) show physical and mathematical schematic diagrams of a 
multi-layer system, respectively. The boundary and interface conditions for each 
layer can be written as follows: 
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Fig. 2.3: A multilayer system (Physical domain) 


> On the boundary (external surface) of the multi-layer system, the heat 
flux or temperature should be prescribed: 


to 

II 

p 

(2.26a) 

II 

(2.26b) 


where Q is heat flux on the surface, Q p is applied heat flux, T is temperature on 

the boundary and T p is prescribed temperature. 

> On every interface Fi between layers j and j+1 heat flux and temperature 
satisfy the continuity condition: 

&=! 

7 ) =1 

where the subscript I indicates I th interface. 
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Layer 1 

' Layer 2 ' 

Layer 3 
Layer 4 



It is important to note that only the perfectly bonded part satisfies the 
above equilibrium conditions. The non-perfectly bonded part, say, an interface 
crack, will be considered as an additional boundary of each layer and the heat 
flux and temperature values associated with them will not be coupled together in 
the final equation. 

The conditions at the interface, I between layers j and j+1 as given in Eq. 
2.27 can be implemented in boundary element equations as follows. 


For the j* layer, the temperature and heat flux are related by. 


T J 


[H‘ Hi] =[G< G/] 

l I _ 


& 

Qt 


where T/ and Qj are the interface temperatures and heat flux of layer j on the 


interface Pj, T‘ and Q‘ the temperature and heat flux on the remaining surfaces. 


Similarly, for the layer j+1, the relationship can be written as, 


r pj + 1 

l H/ * H r l ] r/+ , =[<^' or] 


Q'*' 
Qi *' 


(2.29) 


where T/ +1 and <2/ +1 are the interface temperatures and heat flux of layer j+1 on 
the interface Ti, T J+l and Q J+] the temperature and heat flux on the remaining 
surfaces. 

According to the equilibrium condition (2.27) at an interface, Eqs. (2.28) 
and (2.29) can be rewritten as [6]: 


[> J Hj -Gj] 


T > 

T/ 

Qi 


=[G'][e'] 


and 


H{* G/ +l ] 


jij+l 

T/*' 

Qi *' 


=[G y+1 ][e /+ '] 


respectively. Thus, coupling of the above two equations yields the equation, 

V 


~ H‘ H‘ t 

-g; 

0 

t; 


G‘ 0 1 

~ Q r 

1 

o 

t 

Gj* 1 

H 1 *' 

Qi 


+ 

o 

J 

qM 


jij+i 


(2.30) 


Similar equations can be derived for other layers and the substrate. The 
system needs to be reordered to make the matrix equation in banded form. 

2.6 Application of Boundary and Initial Conditions: 


With the boundary conditions (2.2) to (2.4), the equation relative to the 
boundary points (2.19) can be decomposed as follows: 


(2.31) 


'L h \ j ,f.f T j.f + 1 H2 w T ]F +X(ff2 w -Gl^a^Tjj 


m 


m 


7(3) 


S ^^ij\F,F % j,F + ^ ^^iJ,F,F Q j,F + X1 + 


70) 7(2) 7(3) 

where j(k) means element j is in part r\ of the boundary. Writing the above 
equation for all boundary points produces the system 


H2*f, f T f = G2 ff Q f +S*f 


(2.32) 


For linear time variation, if j 6 Ti or T 2 , 

H2*ij,F, F = H2 i JFF 

S*uf=S uf 

andifjGT 3 , 

Hi i,j,F,F — H 2jj FF G2j j F F OL j 

sV=s IiF +Sffl w ^ 

Consider that the initial temperature field To satisfies the equation 

V 2 r 0 +6 0 =0 (2.33) 

in which b 0 is the value of b at time to. Using the above equation in Eq. (2.1), it 
can be obtained 

(2.34) 

with the initial condition T - To = 0 at time to. 

The boundary integral equation equivalent to Eq. (2.34) is 

c*, (r Ui , r - T Ufi y)[k(T-T 0 )qdrdt = 

(2.35) 

') [k{q-q a )T'drdt+\ [k{b-b,)rdQdt 

<b f 0 

The temperature profile over the entire domain can be calculated after 
knowing variation of temperature over entire boundary using Eq. (2.5). However, 
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this does not require time dependent fundamental solution. A simpler solution 
may be used as given in Appendix C B’. 

2.7 Closure: 

In this Chapter, the basic theory for the transient temperature analysis using time 
dependent fundamental scheme has been presented. Application of the boundary 
and initial conditions has also been given. The modification in the modeling for 
the multilayer systems has also been described. 



Chapter 3 

THERMAL STRESS ANALYSIS 

In this Chapter the temperature field as determined in Chapter 2 will be used to 
investigate a thermo-elastic stresses and failure mechanisms. This formulation 
treats the thermal effects as body forces over the solution domain. The body 
force term is a domain integral which would normally require the discretization 
of the domain into cells. Therefore, these domain integral terms are transformed 
into boundary integrals, preserving the important advantage of the BE technique. 

Special emphasis is needed on the accurate evaluation of boundary 
integrals for sub-micro scale domain. If the materials yield at any point, then the 
stress relaxation by plastic deformation is also taken into consideration. 

3.1 Governing Equation; 

For the case of a material being elastic, homogeneous and isotropic continuum, 
the equilibrium of forces and conservation of energy together with the 
constitutive laws yield the following system of governing equations within the 
linear theory of thermo-elasticity [17], 
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(3.1) 


Mu im +(z + M)u kJd +X, = pd 2 u/dt 2 + r <p. 

(3 ' 2) 

k K c 

in which Ui(x,t) and <3>i(x,t) are the displacement and temperature fields, 
respectively, and %, p. are material properties as given below. For the thermo- 
elastic analysis ® is the deviation of the current temperature T from the value T 0 
corresponding to the undeformed state. The symbols Xi(x,t) and Q(x,t) denote 
the intensities of the body force vector and the heat sources. The relationships 
between various material parameters employed in Eq. (3.1) and (3.2) are: 


Table 3.1: Relationships of material properties 


Parameter 

X 

y 

8 

1/k 

3-d problems 

Plain strain 

X 

Yo 

£o 

1/Kq 

Plane stress 

2m v~ 

1 — o 

„ 1 + L> 

2 LUX 

1 -u 

1 — 2l> 

e ° i 

1-y 

1 1 + L> 

— + e n a 

K 0 1-V 


with 

K = O-M + 3^K £„ = k„= — 

4, P C ,: 

p,, X, -- Lame constant (under isotropic condition) 

a — Coefficient of linear thermal expansion 

Xo — Coefficient of heat conduction 

p — Mass density 

c e ~ Specific heat (under constant strain) 

k — Thermal diffusivity 

u — Poission’s ratio 
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The tensor of elastic constants for a homogeneous and isotropic 
continuum is given by 

= x8„S a + rfSJj, + S„S Jk ) (33) 

The constitutive law that includes deformations and temperature (known 
as the Duhamel-Neumann relation) is given as, 

v,j = 2 M,j + (X£u - Yf)S,j ( 3 . 4 ) 

This can be rewritten as 

Gy ~ C ij]d £ kl ~ Y^ij 

In general, the boundary conditions for the thermal stress analysis can be 
given as: 

> Domain should have prescribed tractions at some part of the boundary, 
i.e. 

T =T onTi (3.6) 

n n 

> Remaining part of the boundary should be constrained, i.e. 

u =u„ onT 2 (3.7) 

3. 2 Boundary Integral Equation: 

The above governing equations can be recast as a three-dimensional boundary 
integral form having body force volume integral (considering no heat source) as 

[ 16 ]: 

«,(P)+ \T,{p,Qyu,(Q)dS(Q)= (3.8) 

S 

\u,j (P, Q) 1 1 {Q)dS (Q) + jU^giXjWdViq) 

s v 

where S is the surface and V is the volume. Ty(p,Q) and Uj,(p,Q) are called 
traction and displacement kernels and can be given for two-dimensional 
problems as: 



(3-9) 


U,j{p,Q) = 


8^w(1-l>) 


(3-4o)ln 


1 


L r(p,Q)_ 


S. | dr (P*Q) friP’Q) 


dx . 


dxj 


t (p n\- 7_I farQ’.fil'lL , w ,arfe0%2) 

' A '’ 4*(l-o)r(p,Q)[ Bn Jf 2 °^ + 2 &, &, 


(3-10) 


+ 


\ — 2v 


4x(l-u)r(p,Q) 


dr(.P,Q) n or(.P,Q) n 


dxj ' dx j 

The body force vector, Xj(q) can be expressed as the gradient of a 
potential function, O(q), as: 

Y( , -aEdtiq) (3.11) 

yW 1-2 v dxj 

Eq. (3.8) is valid for both 2-D plane strain and 3-D generalized plane 
strain problems. The effective coefficient of thermal expansion for two- 
dimensional problems is therefore given by 

a* = a (Plane strain) 

l +u (Plane stress) (3.12) 


a* = a- 


l+2v 


3.3 Two-Dimensional Thermo-Elastic Kernels : 


From the boundary integral equation (3.8), using divergence theorem the body 
force term can be modelled in the following form [17]: 

B, = jM,(p,Q)^Q)<ir(Q) + jN,(p,Q)^-dr(Q) (3 ' B) 

and, the two-dimensional thermo-elastic kernels can be written as follows: 


N ^£^U r2{p ’ Q)Xnr ^ 


(3.14) 

(3.15) 





These kernels can be further expanded to 


(3-16) 

(3.17) 


3. 4 Numerical Implementation: 


The evaluation of the kernels is achieved by dividing the surface (boundary) into 
elements and performing the numerical integration over each element. Therefore, 
the discretized form of the BIE can be written as : 

ML I ML 


m=\ c=l _i m=l c=l _1 


m=l <?=1 


-1 


m-l c=I 


(3.18) 


where M is the total number of elements and L is the number of nodes per 
element (eight for three-dimensional quadratic elements). Taking each node in 
turn as the load point P and performing the integrations, a set of linear algebraic 
equations emerges as follows: 

(3.19) 


[4“H*]['M^Mg][£' 


where the matrices A and B contains the integrals of kernels Uy and Tij, 
respectively, while the matrices F and G contain the integrals of the thermo- 
elastic kernels, Mi and N[, respectively. 

As mentioned earlier, the temperature and temperature gradient 
calculated from transient analysis must be substituted in Eq. (3.28). No further 
manipulations are necessary to form the solution matrix, because the thermo- 
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elastic contributions are simply added to the right hand side of the solution 
matrix, as follows: 

MM=MM+[/M*] (3 - 20) 

where [x] contains all the unknown functions (whether displacements or 
tractions), and [y] contains all prescribed values (whether displacements or 
tractions). 

3.5 Boundary and Interface Conditions: 


Assembly of the [A], [B], [F] and [G] for the multilayer can be done in the 
similar way as descried in the section 2.5. The boundary and the interface 
conditions now can be given as: 

> On the boundary (external surface) of the multi-layer system, the traction 
must be given, such as: 


'„=U=? (3.21) 

where, tn and t t are the normal and tangential components of the traction, T n 
and T t the applied loads in the normal and directions, respectively. In 
addition, the multi-layer system should be constrained (with specified 


displacement) at some other location of boundary. 

> On every interface Fi between layers j and j+1 normal and tangential 
traction and displacements from both sides of the interface satisfy: 

t, 


In 

U 


= / i _ t j + 1 

l In l ln 


t J ' =-t j+1 
l It l It 


~ uj = u J+l 


... y_ y+i (3*22) 

’’In — “In M '/n ’ U p ~ U It ~ U It 

where the subscript In indicates normal traction component at interface (I), and 
subscript It indicates tangential traction component at interface (I). 
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External Surface 





Fig. 3.2: A multilayer system (physical domain) 

An advantage for the BEM as applied to multi-domain problems is that 
both traction equilibrium and displacement compatibility conditions between 
layers are explicitly satisfied, as compared to FEM, in which only the 
displacement compatibility condition is explicitly satisfied [13]. 

3.6 Modification for Nearly Singular Integral for Sub-micro scale 
domains: 

The difficulty in applying the BEM for thin-structure and crack problems is the 
evaluation of nearly singular integrals. The integrals in BEM, which determine 
the influences matrices, contain singular kernels of the order 0(1 /r) and 0(ln r) 
in 2-D elasticity case, where r is the distance between the source point and the 
integration point on the boundary element. When the source point is very close 
to, but not on the element (0 < e « 1, where e being the distance between 
source point and element of integration), although the kernels are regular in 
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mathematical sense, the values of the kernels change steeply in the 
neighborhood of the source point [12]. The standard Gaussian-quadrature is no 
longer practical in this case since a large number of integration points are 
needed in order to achieve a required accuracy. 

As the kernel functions are given in Eqs. (3.9), (3.10), (3.15) and (3.16), 
the displacement kernel Tjj contains (1/r) term and ln(r) terms, which are 
singular when the distance r approaches zero. Because of existence of these two 
singular terms, the integrals on a typical element T, 

(3 ' 23) 

r 


and. 


ju,( P> Q) tj (Q)dr(Q) 


are called weakly (the Logarithm type) and strongly singular integrals, 
respectively. In the modeling of thin structures in the micro or nanoscale by 
BEM, one possible way is to evaluate them by singularity subtraction and 
nonlinear coordinate transformation [14]. 


3. 6. 1 Regularise nearlv-singular integral by singularity subtraction: 


Rewriting Eq. (3.23) in the following form: 


r 

+ 


(3.25) 


T(0 


where, Q 0 is the closest point on T to p Fig. 3.3. As Q — ► Q 0 , the term 
[u .{Q)-u {Q )] has the order of °( r ’) where r = a§- Then the order of the first 

integral of right hand side in Eq. (3.35) reduced to 0(r’)/0(r). This integral is 


nearly weakly singular when e is very small. The last integral on right hand side 
of Eq. (3.25) can be evaluated using the function value at the two end points as: 

fo C p, Q)Uj (Q) dT(Q): = C,C 2 ( e v e 2j - e 2l e tJ )]n, rP 
r 1 (3.26) 

+C, {C 2 S : fi + s u s u (0 + hm20)~(S u S 2j + )cos 26 

1 

+S 2i S 2 i{0--s m20)} 

2 <h 



Fig. 3.3: Source point near to the line T of integration 
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p 

Fig. 3.4: Element T in natural coordinate system 


Now, Iet<£ ~r/ m (m>2), than the above equation reduces to: 


2 1 



tj +e 


drj 



rj +e 


drj 


(3.28) 

I 


where m is the order of nonlinear coordinate transformation. When m is large 
enough, the first integral of right hand side of Eq. (3.25) can be evaluated 
accurately. For quadratic elements, when e is on 10^ scale or smaller, m = 6 is 
accurate enough. But best result can be obtained when m lies between 12 to 14 
using 10-20 gauss points. 


3.6.2 Evaluation of nearly weakly singular integral: 

In 2-D elasticity problems, the weakly singular integrals given in (3.24) contain 
the singular term ln(r). In the singular case, if r = 0, a special logarithmic 
Gaussian quadrature can be employed. However, when r is very small but not 
zero, this nearly weakly-singular integral can cause difficulties in the 2-D BEM 
procedure. Similar nonlinear coordinate transformation as developed in the 
previous section is used to evaluate the nearly weakly-singular integrals 
accurately. 
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Considering the natural coordinate system again, with nonlinear 
transformation,^ = rj m , the integration of the term in Eq. (3.24) can be written as: 

\\nrdr = 2 J-ln(£ 2 + e 2 )d<* ( 3 - 29 ) 

r o ^ 

i 

= j j mrj m ~ 1 \n(Tj 2m +e 2 )drj 

o 

The above integral can be evaluated accurately using a small number of 
Gauss points when m is large enough. 

3. 7 Stress Relaxation: 

The residual stresses at any point inside aluminium layer derived from the 
previous Sections may exceed the yield limit. Then, these residual stresses will 
be relaxed either by plastic deformation or by diffusion of voids. The stresses, 
thus, redistribute themselves according to the temperature profile and material 
properties locally. 

The relaxation of stresses due to diffusion of voids can be taken into 
account by the use of addition eigenstrain (principle strain) term Sij R in the 
thermal strain [8]. The total eigenstrain sy** can be given by, 

£~=e T +e« (3.30) 

4=* r +4 (3.31) 

4=^+4 = * r -(4 +4) ( 3 *32) 

where £ T denotes thermal strain. 

Eq. (3.32) is obtained from the condition that the volume of the 
aluminium line does not change during the relaxation process: 



The stress after relaxation can be calculated from eigenstrains as 


follows: 


2+r 


2*o 


or n — ( r-£ n ^ — H £,,) 

1-0 (1+r) 2 11 (1+r) 2 22 1 +r 33 ' 

r(2r + l) ** 2 *r*v 


, r 

^22 ~ ~ 7 ~~ ( + 


“f”* 


l-u\l + r) 2 11 ' (1+r) 2 ~ 22 


1 + r 


cr 33 — 


_ 2 n 


( u (, , £ U + 1 . £ 22 ) £ 33 ) 


l-o 1+r 


1+r 


(334) 
4 ) ( 335 ) 

(336) 


where r is the aspect ratio of aluminium layer. 

Unless the constraint by the surrounding matrix is removed the stresses 
in the aluminium line can not become zero even after complete relaxation. These 
stresses can be found in the following manner. The elastic strain energy per unit 
volume E 0 can be expressed as: 

(3.37) 

where the usual summation convention over repeated indices is adopted. From a 
thermodynamic point of view, the change in the elastic strain energy under a 
constant temperature is the change in die Helmholtz free energy of the material. 
Therefore, the relaxation process will proceed until Eo takes it’s minimum value. 
The stresses after complete relaxation can be obtained by applying the energy 
minimum conditions, i.e. 

•?4 = 0and-^- = 0 (3.38) 

5s" de» 

The equilibrium values of e*, corresponding to the complete 

relaxation can be obtained as 

T _ ( 1 + *^)[(4r 2 - r - 2) - 2(2r + l)(r - 1 )o] 

€u (2r 2 + r + 2) - 4rv — 2(1 - r 2 )t> 2 

~T _ (1 + t>)[-(2r 2 - r - 2) + 2(r + 2)(r - l)t>] ~ 4Q . 

22 (2r 2 +r+2)-4ru-2(l-r> 2 
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(3.41) 


"33 


and, the stresses will be 

_/2 , R R 

CTi j — <7 22 ~ CT 33 — CT — — 6/2 


2(1 + *j)[(r 2 - r + 1) - (r - 1) 2 u j 
(2r 2 + r + 2) - Arv — 2(1 - r 2 )v 2 

r{ l + u) 


(3.42) 


(2r 2 + r + 2) - 4ri> - 2(1 — r 2 )u 2 
It is interesting to note that internal stress becomes hydrostatic after the 
complete relaxation. 


3. 7. 1 Stresses after plastic deformation: 

Relaxation by plastic deformation can occur instantaneously by the conservative 
motion of dislocations. Von Mises or Tresca criteria can be applied to predict 
the onset of the plastic deformation under the condition of triaxial stresses. For 
the simplicity of the Tresca criterion, this criterion is being applied here to 
discuss the plastic relaxation of the aluminium line for 0 < r < 1 . According to 
the Tresca criterion, yielding occurs when 

a„-<T 21 >cr y (3.43) 

where er 33 is largest and cr 22 is smallest among all the stresses and cr y is the 

uniaxial yield stress. The stresses after the plastic relaxation can be evaluated by 
applying the energy minimum conditions under Tresca criterion. From Eq. (3.43) 
we have 

° P ii- G y ( 344) 

The superscript P indicated values after plastic deformation. Similarly 
the relation between s p 12 and s F u can be obtained by using Eqs. (3.34) to (3.36) 
and (3.44), following similar procedure as described for diffusion relaxation, as 


e p 22 — Cj +C 2 e ,, 


(3.45) 
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where 


; = S+/) 2 ( ^v) r , i-i> , 

4r 2 (1 - u) + r(5 - 4l>) + 2 ^ 1+r // y 

C 2(l + r) 2 (l-t>) + r 

2 4r 2 (1 - l>) + r(5 - 4u) + 2 




From the minimum condition, —— = 0 , e xx can be obtained at the end of 


de ; 


the relaxation by plastic deformation as: 


[C 2 (2C 3 +C 2 C 5 ) + C 4 ]^ +(C 6 +C 2 C 7 )£ r +C I (C 3 +C 2 C 5 ) = 0 
where C 3 to C 7 are 

C _ 2r 2 (1 - 1 >) + r (5 - 4u) + 2(1 - u) 

3 (1+r) 2 

r - 2r 2 (l-L>)+r(5-4o)+4(l-u) 

(1+r) 2 

C 2r(l+L>) 

6 0+r) 


(3.46) 


_. 2q+o) 

7 (1+r) 

is calculated from Eq. (3.46). ^ can be calculated from Eq. (3.45) with e xl . 
can be obtained from 

4 =-( 4 + 4 ) (3.47) 

The stress field after the plastic relaxation can be solved by putting 
above calculated eigenstrain into equations (3.34) to (3.36). 

In this section it has been found out that the internal stress a* can, in 

principle, be reduced to a R if the complete relaxation is possible. As the 

relaxation process can occur either by the plastic deformation of aluminium or 
by the diffusional atom movement of aluminium. However, at low temperatures 
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where the diffusional relaxation is practically inoperative, only the relaxation by 
plastic deformation can take place [6]. So, the stress relaxation by complete 
relaxation can be neglected at low temperatures since it will be of very low order. 

3. 8 Closure: 

In this Chapter, the thermal stress formulation has been explained briefly. 
Modification in the kernels due to the thin multilayer structure application has 
also been presented. Stress relaxation by plastic deformation has also been 
modelled. 



Chapter 4 

THERMAL GROOVING 


Recent works have revealed that stress-induced migration causes serious 
damage to the aluminium conductor. The failure caused by thermal stresses is 
classified into two types depending upon the temperature. (1) Low-temperature 
failure: (at temperature lower than 273 K) the cyclic stress leads to fatigue 
fracture. (2) High-temperature failure: (at temperature higher than 373 K) the 
stress causes motion of aluminium atoms and the diffusion causes micro defects. 

In this chapter, a stress-induced migration, i.e. diffusion along the 
surface and grain boundary is being modelled. In addition, the model is further 
modified to define a fracture criterion so that it can be useful for prediction of 
time to failure. The effect of the stress concentration near the crack tip is also 
modelled using empirical relations. 

4.1 Introduction: 

Aluminium when deposited on a surface followed by cooling forms a 
polycrystalline structure; the size of the grains thus formed depends on the 
growth conditions. Fig. 4.1a shows an ideal section of such a grain boundary- 
surface intersection. The exact shape of the groove is determined by the balance 
of local thermodynamic forces (the two surface tensions and the grain boundary 
tension) requiring an equilibrium angle of qr, as it shown in Fig. 4.1b. With the 
development of the groove, ridges are also formed on the surface lying on either 
sides of the interface. These ridges try to flatten by surface diffusion; however 
this flattering process affects the equilibrium angle and makes the groove to 
deepen (Fig. 4.1c). Here, this situation is prior to deposition of the passivation 



coating so that there are no residual stress along aluminium layer, thus as a 
result there is no atomic flux along the grain boundary. The groove developed in 
this way will initially grow very fast and then slows down to reach an 
equilibrium situation in a period of few hours. This process is known as thermal 
grooving phenomenon [8]. 

Thermal grooving alone, therefore, does not cause the failure in the 
aluminium line by itself, but it provides an early nucleus for grain boundary 
cracking when wire is subjected to high stress levels. The stress perpendicular to 
the grain boundary activates the atomic flux along the boundary which enhances 
the grooving and ultimately resulting in open-circuit failure. 

4. 2 Estimation of Thermal Grooving: 

A typical aluminium conductor in a microelectronic package is shown in Fig 
4.2a. In this structure the width of the aluminium layer is so thin that the wire 
can be considered as a sequence of grains defined by a periodic intergranular 
boundary which completely traverses the width of the wire. Although the 
bamboo like structure is effective to prevent electromigration, stress-induced 
migration inflicts serious damage on the conductor. In the present work a typical 
layer width of 1 pm is considered which is under a remote tensile stress. 

Using the symmetry of the problem, only one quadrant of the groove is 
analyzed. The initial surface is discretized by a set of N equally spaced points as 
shown in Fig. 4.2(b). 

As explained in reference [4, 7] the key relations to analyze the thermal 
grooving phenomenon are as follows: 


38 



(a) Semi-infinite body with flat surface 
and grain boundary 


Grain boundary 

Surface 


(b) Thermal grooving due to 

Radius of curvature - 1* surface diffusion 


(c) Grain boundary grooving due to combination of 
diffusion along the surface and that along the grain 
boundary 


Fig. 4.1: Explanation of grooving due to diffusion along the 
surface and grain boundary 






2W = 1 jim 


Grain 
B oundary^\| 


r« — H 


2L = 20 jum 


(a) Model conductor with 
bamboo like structure 


Z 



(b) Discritized quadrant of 
the conductor surface 


Fig. 4.2: Modal conductor 
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The normal velocity of the atoms (Fig. 4.3) at the surface relative to the 
adjoining solid material can be written as: 

v„ = ~(D s S s Qr s / kT)d 2 tcJ ds 2 (4.1) 

where, D s is the surface diffusivity, 8 S is the thickness of diffusion layer, Cl is the 
atomic volume, y s is the surface free energy, k is the Boltzmann constant, T is 
the temperature, k is the curvature of the surface, and s is the arc length (Fig. 
4.4). 

From the equilibrium of the local thermodynamic forces, angle y, as 
shown in Fig. 4.1b, at the interface of the grain boundary and the free surfaces, 
can be given by: 

2/s cos (p - y b (4.2) 

where, y b is the free grain boundary energy. The curvature at the groove 
tip, K t ip is defined as [7]: 

** ^KpJr.-Vm-atWiFidK/as^ (4.3) 

where Ki is the stress intensity factor in mode-1, cr*, is the far field stress, a is 

the groove depth, W is half width of aluminium layer, and F is diffusion ratio. 
Stress intensity factor for double-edge-cracked-semi infinite plate under uniform 
tension can be given as [21]: 

K, =crJxa-F l (a), a = ^ - (4 15) 

(a) * 1 . 12 - 0.23 la + 10.55a 2 - 21.72a 3 + 30.39a 4 

, and, 

F = D s S s l D h 8 b (4.16) 

The value of F usually depends not only on the material but also on the 
temperature and the microstructure of grain boundary and surface. Here, D b is 
the grain boundary diffusivity; 8 b is the thickness of diffusion layer at the grain 
boundary. 
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Fig. 4.3: Groove surface profile at time f and f+1 



Fig. 4.4: A typical node i subjected to two adjacent elements 

Curvature at the other the other nodes can be calculated by the technique 
given in Appendix ‘D\ Using the curvature the velocity at each node can be 
calculated except the groove tip. The new profile of the groove surface can be 
calculated by using the velocity at each node. The new coordinate of the groove 
tip can be obtained by, 
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(4.17) 


V=o 

The new and old position of the groove tip can be used to calculate the 
velocity of the groove tip v N as follows (Fig. 4.3): 


N 


At 


(4.18) 


where At is the time step used. The implementation details for the thermal 
grooving are given in Appendix ‘D’. 

Some of the assumptions made while deriving the above relations are: 

1 . Creep cavity growth is uniform along the grain boundary. 

2. The properties of the interface are independent of its orientation with 
respect to the adjacent crystals. 

3. The only mechanism operative in the transport of the material is surface 
and grain boundary diffusion. 

4. Effects of the inplane stresses are negligible in thermal grooving process. 


4,3 Closure: 


In this chapter, stress induced failure of aluminium line at high temperature in 
terms of surface and grain boundary diffusion has been modelled. To find out 
stress at the crack tip use of mode-1 stress intensity factor has been made. 
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Chapter 5 

RESULTS AND DISCUSSION 

The present BEM method has been applied to perform failure analysis of an A1 
line of multilevel VLSI circuits. Two ways of failure (plastic flow and thermal 
grooving) were considered. The failure analysis of the A1 interconnectors both 
during the processing and after the processing is presented here. The results are 
organized as follows: 

1. Failure analysis during processing 

2. Thermal grooving 

5. 1 Validation of the BEM Code: 

To validate the BEM code an example problem has been done. In this problem a 
square plate of silicon with lmxlm dimensions was considered as shown in Fig. 
5.1. A heat flux of 1000 W/m was applied on one of its face. Other faces were 
left free to dissipate heat by natural convection. The problem was analyzed 
using the same BEM model discussed in previous Chapters. The transient 
analysis has done till the structure reached the steady state. Property of the 
silicon is given in Table 5.1. The results have been validated through finite 
element software package ANSYS. 

The temperature variation at the middle node of the different surface 1, 
2 and 3 as shown in Fig. 5.1 is plotted with time in Fig. 5.2. The differences in 
the three plots are due to the time taken by the heat flux to reach at that location. 
This plot indicates that initially the temperature increases faster and slows down 
later. Comparison of the temperature variation with time obtained from the 
presented BEM technique and finite element software ANSYS at middle point 
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of surface 1 is shown in Fig. 5.3. The results obtained by the BEM code differs 
by 1.5% to the results obtained from ANSYS. 

Fig. 5.4 shows the temperature profile over the entire domain at steady 
state. This shows that the temperature is highest at the middle node of surface 1 
and lowest at surface 3. It shows that temperature is maximum at the point of 
application of heat flux and reduces further away from it, which is true 
qualitatively. Fig. 5,5 shows the temperature variation over entire domain 
obtained from ANSYS software. The trend of temperature variation over the 
domain by ANSYS matches with the trend of temperature variation over the 
domain by BEM code. 

Fig. 5.6 shows the total heat flux over entire domain at steady state. Total 
heat flux is maximum towards the comer of the surface 1 and minimum is at the 
middle region of surface 3. 

5.2 Failure Analysis of Aluminium Line d urins Processing: 

The present analysis is done for a 4-layer printed circuit board (Fig. 5.7). The 
fabrication process in this work is accomplished in two steps. First, the silicon 
buffer is covered by borophosphosilicate glass (BPSG) at 673°C, and cooled to 
room temperature. In the second step, the 1pm width A1 line is formed over the 
BPSG, and passivated with phosphosilicate glass (PSG). Deposition temperature 
of the PSG passivation layer was also 673°C. The complete structure is cooled 
down to the room temperature. Properties for the above four materials are listed 

in table 5.1. 

The problem is treated in a semi coupled manner: for each time step a 
thermal analysis is performed, then the resulting temperature profile over the 
boundary are used for the stress field estimation. Only one half of the domain is 
analyzed due to the symmetry of the problem. The geometry and the boundary 


45 



conditions of a typical multi-level structure are shown in Fig. 5.7. The cooling 
process is controlled by natural convective boundaries with appropriate values 
of convective boundaries coefficients. The ambient temperature is assumed to be 
25°C. A typical time step of 0.0001 second is used in the present work. It is 
assumed that only aluminium has plastic behavior and silicon and glasses are 
treated as elastic materials. To adequately simulate the large size of the wafer, 
the silicon wafer used was assumed to be 10 times the thickness of the BPSG, A1 
and PSG combined. All nodes falling on the line of symmetric are constrained. 


Table 5.1: Property data for different layers 



Silicon 

BPSG 

Aluminium 

PSG 

E (GPa) 

130 

120 


79 

p (Kg/m J ) 

3330 

2340 

2700 

2300 



u 

0.28 

0.2 - 

0.33 


a 

2.3e-06 

l.le-06 

23.5e-06 

l.le-06 

Cp (J/Kg*K) 

712.25 

787.5371 

904.322 

756.471 

h (J/m z *K) 

21 

4 

56 

7 

X(J/m*K) 

150 

10.4 

235 

12 


As explained earlier, the structure is analyzed in two steps: 

Step I: In this step BPSG layer is deposited over the silicon layer at 673°C and 
structure is cool down to the room temperature. The discretized domain for step- 
1 is shown in Fig. 5.8. The cooling curve for the first step is shown in Fig. 5.9. 
This curve shows that initially the cooling rate is very fast and decreases further 
as the difference between the layer temperature and atmospheric temperature 
decreases. It also indicates that the silicon and BPSG layer takes only two 
seconds to cool down. It was also observed that the temperature profile over 
entire domain at any instant is quite uniform. This is due to the fact that the 

























structure is of micro-scale and there is not much difference in temperature of 
boundary and interior points. 

The residual stress CT ^ , in the silicon buffer and BPSG layer at the end of 

Step-1 is shown in Fig. 5.10. It is clear from the stress profile that the silicon 
layer is under tension while the stress in the BPSG layer is compressive. This is 
due to the fact that the thermal expansion coefficient of the silicon is almost two 
times that of BPSG. The stress is varying parabolically along the interface 
having its peak at the middle of interface. 

The residual stress <j w plot is shown in the Fig. 5.11. It can be seen from 

the stress profile that the maximum stress is less than that of the cr M . Further, 

the a is of tensile nature in the middle region of the interface, while it is of 

yy 

compressive nature towards the either ends of the interface. The residual stress 
r xy , plot after step 1 is shown in Fig. 5.12. The stress profile indicates that the 

r at the interface is of tensile nature towards the free end of the interface, 

while it is of compressive nature towards the line of symmetry. 

The residual stress cr a plot is shown in Fig 5.13. Comparing the stress 
distribution of and cr a , it can be viewed easily that the stress values cr a is 
more. 

Step-2: In this step aluminium and PSG layer is formed over the BPSG layer at 
passivation temperature of 673°C and full structure is cool down to the room 
temperature. For Step-2 discretized domain is shown in Fig 5.14. The stress 
profile calculated from Step-1 is taken as initial stresses for this Step. Fig 5.15 
shows the variation of temperature with time at the middle of the. 

The residual stress after deposition of PSG passivation layer and 

plastic relaxation is shown in Fig 5.16. The stress field obtained shows that the 
aluminium layer is under tension and the maximum stress lies in the aluminium 



line towards the PSG layer. This is due to very high value of thermal expansion 
coefficient of aluminium with respect to other layers. 

Similarly, the residual stress and <j zz plots after plastic relaxation at 

the end of step-2 are shown if Fig 5.17 and 5.18, respectively, cr^ is maximum 

in the aluminium layer while it is of almost zero magnitude at silicon layer. 
Although there is no thermal and structural analysis performed in the axial 
direction, the plain stress model used can predict the stress field in this direction 
also. As shown in Fig. 5.14, the cr^ has values closer to that of cr^ . 

5.3 Thermal Grooving: 

According to Eq. (53) groove growth depends on the axial stresses appearing in 
the aluminium line. Fig 5.19 shows the charge carrier part of the aluminium line 
subjected to the cr a . This plot reveals that the stresses in the axial direction do 
not very much, so it can be assumed that the far field stress to be constant 
throughout the aluminium layer. The value of the far field stress used in present 
work is the average of stresses appearing in aluminium layer. The values of the 
parameters required to analyze thermal grooving process is listed in Table 5.2. A 
typical time step of one second is used in the present work. 

The modelled groove profile is shown in Fig 5.20. Fig. 5.21 shows that the 
initial surface profile and the final surface profile at the grain. The final profile 
here is related with situation when the half width of the aluminium layer is 
reached by groove tip and this corresponds to the open circle failure of the 

aluminium layer. 



Table 5.2: Property data required to analyze thermal grooving process 


Ys 

1.5 J/m 2 

D b 8 b 

2.64e-23 nvVs 

Q 

1.66e-29 m J 

k 

1.38e-23 J/K 

F 

100 

Yb/2 Ys 

1.6 


Fig 5.22 is plotted between groove depth and time. Initially groove tip 
moves faster but slows down as the crack move forward. Fig. 5.23 shows the 
varation of non-dimensionalized groove depth with nondimensionalized time. 
The factors used to non-dimensionalize the width, time and stress are given as 
follows (Section 4.3): 


T = (D s SprJ)/(kTW > ) 

(5.1) 

s=s/W 

(5.2) 

x-xlW 

(5.3) 

a -a/W, 

(5.4) 

v n =v n km 1 KD t s s nr,l 

(5.5) 


where values for all the property data is given in Table 5.2. 

Fig 5.24 is plotted between groove growth rates (i.e. velocity at the crack 
tip) and time. The velocity decrease at the early stages while it becomes constant 
afterwards. The decrease in growth rate at early stages is governed by the 
process of thermal grooving due to surface diffusion. Fig 5.25 is plotted between 
non-dimensional groove growth rates (i.e. non-dimensiona I velocity at the crack 
tip) and non-dimensional time. 

Fig 5.26 shows that how stress intensity factor is varying with time. It is 
initially increases faster and it slows down as the crack moves forward. 
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The time to failure for the aluminium layer as calculated by the analytical 
model are few seconds [5] while the experiments results have shown that the 
value of time to failure is typically about hundred hours. The presented 
numerical method also predicts that the time to failure of aluminium layer is 
about hundred hours (Fig. 5.27). This is due to the fact that analytical model 
assumes that grain boundary diffusion process is much faster than surface 
diffusion (as shown in Fig 5.28), and, under the application of far field stress, 
grain boundary diffusion control the time to failure. However, the presented 
numerical method predicts that two processes exist simultaneously, which is 
more realistic. The two processes are coupled at the groove tip by continuity of 
chemical potential and atomic flux. 
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Fig. 5.1: Test Problem 
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Fig. 5.3: Comparison of the Temperature variation with time by developed BEM 

code and ANSYS 


Fig. 5.4: Temperature profile over the domain after reaching the steady state 
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Fig. 5.9: Variation of temperature at mid point of the interface with time 
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Fig. 5.10: cr w stress field distribution at the end of first step 
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Fig. 5. 11 : cr stress field distribution at the end of first step 


z stress field distribution at the end of first step 
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Fig. 5.14: Discritized domain for second si 
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Fig. 5.15: Temperature variation at mid point of aluminium layer with time 


stress field distribution in the final structure after step 2 
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Fig. 5.17: a stress field distribution in the final structure after step 2 


stress field distribution in the final structure after step 2 
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Fig, 5.23: Normalized groove depth Vs normalize time 
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Fig. 5.25: Normalized groove growth rate Vs normalize time 



Fig. 5.26: Stress intensity factor Vs normalized time 
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Chapter 6 

CONCLUSIONS AND FUTURE SCOPE 

This chapter summarizes the work, and will set forth goals for future work. 

6.1 Conclusions 

In this work, the residual stress behavior in the aluminium interconnections in a 
multilevel printed circuit board during processing has been analyzed. Boundary 
Element method has been used for the purpose with required modifications to 
calculate the nearly singular integrals accurately in the submicron scale domains. 
A simplified model for the stress relaxation by plastic deformation has been 
applied* Stress-migration failure caused by diffusion along the free surface and 
grain boundary has been modelled accurately. Effect of the stress concentration 
near the crack tip has been taken care of. The following are the main 

conclusions drawn from the present work: 

> Use of boundary element method is promising in the thin film 

applications. 

> Firstly, a multilevel structure can be modeled step by step, following the 
manufacturing process. 

> In the initial step of fabrication, laying BPSG layer on silicon wafer, 
is high at the interface of BPSG and silicon layer. 

> In the second step of fabrication, lying passivated aluminum layer on the 
BPSG layer, a xx and are high in the aluminium layer. 

> Wherein stress zones near the interface of aluminum and PSG layer are 
beyond yielding which get relaxed by plastic deformation. 



> The surface diffusion is the important mechanism which causes thermal 
grooving. 

> The groove will propagate through half the thickness of aluminum layer 
in 200 hours causing failure. 

> The used step by step technique can be readily extended to more 
complex situation and material combinations. 

> By the use of time dependent fundamental scheme, computation cost 
gets reduced. 

> The result obtained for the thermal grooving give the qualitative 
agreement with experimental observations. 

6.2 Scope for Future Work: 

The present work provides some insight to the residual stress behavior and 
thermal grooving process in aluminium layer of a multilayer printed circuit 
board. The following work is also of interest and can be pursued as a future 

work: 

> A more realistic model for plastic deformation can be applied. 

> Possibility of nucleation and growth of a crack inside aluminium line can 

be investigated. 

> Effect of the any singularity present in silicon wafer like crack, inclusion 
can be investigated. 

> Failure of aluminium line due to electromigration can be mvestigated. 
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Appendix A 


A. 1 Evaluation of the Coefficients of Matrices HI, H2, Gl, and G2: 

It is important to notice that the fundamental solution T* and its derivative q* are 
function of x; thus, in the linear case, we have that H 1, iP ,-H 1 . . if F -/= 

iyj\F yf J 

F -/, and the same for H2, Gl and G2. This means that these matrices can be 
stored and reused whenever needed. Thus, for computer efficiency, it is possible 
to compute and store all these matrices, which depend only on geometry and 
time step values, and perform several simulations, for instance for different 
boundary conditions. 

For nonlinear problems, such preliminary calculation is not possible 
since the diffusivity k is temperature dependent. On the other hand, one should 
avoid recomputing all coefficients for each new time step due to efficiency 
reasons. This problem can be overcome by the use of 4N 2 functions, namely 
Kly , K2y , K3y , K4jj , closely related to the coefficients Hly^f , 5 

Glij,F,f , and G^j^f but depending only on the product kx. 

A.l.l Time Integration: 

Expressions (2.6) and (2.7) can be written in the condensed form 

r =Cq ' = Wy^ exp ^ (AA) 

with, 

0 = —^ and £- = 0 if X' =T 

IT' 2 
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and s - 1 


if X*=q 


2d 

" ^ /2 

After inverting the order of integration in expressions (2.14) to (2.15), 
the following time integrals have to be evaluated. 

v 

71= \kX'dt 

V-l 

V 

12= J ktX'dt 

V-l 

For the linear case, k is a constant; for the nonlinear one, it is assumed 
that the time step At is sufficiently small for the approximation k = constant 
between tf„i and tf to be valid, e.g. k(t) = k(tf.i) over Atf. 

A. 1.2 Evaluation of II: 


Substituting the expression for X*, i.e. 

k 


n=c | 


s/2+e 


( 4*0 

and performing the change of variables 


exp 


( ^ dt 


with 


dx 


4 kr Ak{t F -t) 


r . x . 
-dt — —dt 


Akx x 
The integral II can be rewritten as 

kx e 


1\ = C J 


V-l 


(4*r) 


s/2+e 


■dx 


(A.2) 


(A-3) 


(A.4) 


(A.5) 


71 



Eliminating the product kx, according to expression 2.20, gives 


71 = 


C 


4 r 


The result of which is 


71 


s+2s-2 


:‘ ,2 *‘- 2 e - x dx 


X S - 1 


4 r 


s+2s—2 


r -+€-\,x 
2 


where T is incomplete Gamma function, 

OO 

T{n,x)=\e n -'s- e de 

X 

A. 13 Evaluation of 12: 


7-i 


(A.6) 


(A.7) 


Substituting the expression for X*, i.e. 


3 (4 hr) 

7~» v 

and taking into account that x = t F - 1, 


I 

12 = c] — exp' 4fa> dt 


I2 = t r n-C | ,„^ -exp ( 4fe> <# 


, (4*r) 

7-i 


=f F 71 — J2 


(A.8) 


(A.9) 


Introducing again the change of variables given by (A.l), the integral J2 
becomes 


J2 = C J 


hr 1 


\s/2+s 


3 (4kr) s x 

7-i v 


-dt 


(A. 10) 


or, eliminating the product kx, 


J 2 


C 1 


16r 


5+2f-4 


J X 


,s/2+s—2 -x 


e X dx 


(A. 11) 


x r - 1 
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This gives 


1 6r s+2E ~ 4 k 

A.1.4 Evaluation ofr(n,x): 

The following recurrence relation can be used for evaluating the incomplete 
Gamma function, 

r(n+l > x)-riT(n,x)+x n e~ x (A.13) 

T(0,x ) = E l (x) 

r(^, x) = Jx]\-erf{4x^ 

r(\ y x) = e~ x 

in which Ei is the exponential-integral and erf the error function. 

A.1.4. 1 Numerical integration of the HI, H2, Gl, and G2_: 

The final step in the calculation of the coefficients of matrices HI, H2, Gl, and 
G2 is the spatial integration along the boundary elements. Calling 

z = kr = k(t F -t) 

It is possible to introduce four functions of z, namely Kl ij ,K2 iJ ,K3 i j , 
and KA i J such that 

J/i( 7 ;-)rfr=[ A, (A.i4) 

r , 

\l\(q’)dT = [K2j'; (A- 15) 

r / 


(s \ 

1- £ — 2,JC 

U J 


(A. 12) 
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(A. 16) 


\j2(T;)dT = Uk3J'; 

p Kr 

j 

\j2(q-)dV = (A. 17) 

p /V 

1 / 

The preliminary calculation of the above 4N 2 functions (N: number of 
boundary nodes), for point distributed between 0 and y max , makes it possible to 
determine by interpolating the values of K\ . to KA tJ at y/, for all / and all 


pairs i, j, during the resolution procedure. 

Thus taking into accounts expressions (2.12) to (2.15), from the 
definition of II, and J2: 


-5-{«/ -«i> 2 >,ir 4^ir} 


(A- 18) 
(A. 19) 
(A.20) 
(A.21) 


The boundary elements employed in the present work are super 
parametric, i.e. the shape function used to model the geometry are of higher 
order than the interpolation functions used to represent the functional variation. 

For two dimensional problems, the expressions for K1 to K4 are 
(omitting the subscripts i, j for simplicity): 


+i 


Ak 


] E \{ x )\A dr > 

-i 

(A.22) 

:]4 e'W* 

-l r 

(A.23) 


1A 



(A.24) 


E3 = — fr 2 E.(x ) \j\drj 

Ujv\ x 'J ' 

K4 = — f dE x (x)\J\dri (A.25) 

8^r_* 

where r| is the natural coordinate and \J\ the Jacobian of the transformation 
from F to q, i.e. 



Calculation of the off-diagonal coefficients of matrices HI, H2, G1 and 
G2 involves only regular integrals, for all time levels tf up to and including the 
actual time tF. In this case, it is noted that x = 0 and x — *■ oo according to (A. 1); 
thus, all integrands in expressions (A.22) to (A.25) become null. 

For the diagonal coefficients (i.e. i = j in expressions (A. 18) to (A.21)) 
we have that, for any t f t F , there is a singularity in the calculation of K1 at the 
source point I where f = 0 (and x = 0). It is possible to remove this logarithmic 
singularity by a change of variables 

r] = sgn(^) 

in which sgn(£) takes the sign of § and the integration limits are not changed. 
The Jacobian of this transformation is, 

It can be seen that this Jacobian is zero at the singular pointy tj = § = 0, 
thus the new integral is regular and can be evaluated using standard quadrature 
schemes. 
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Appendix B 


B.l Fundamental Solution for Steady State Temperature 
Calculations: 


To calculate two dimensional temperature profile over the domain, the 
governing equation is Laplace’s equation, which can be written as: 


2„-. d 2 T d 2 T 


VT 




0 


dx 2 dy 2 

which can be recasted in the boundary integral form as follows: 


c(p)T(p) + \K,(p,Q)T{Q)dT(Q) = \K^p,Q)^^-dT{Q) (B.2) 


where Ki and are the kernels. These kernels can be evaluated as: 


(B.l) 


K,(P>Q) 


3^P,Q) 


8n (B.3) 

k 2 (p,Q)=Hp,Q) 

where X is the fundamental solution. X, for the steady state thermal analysis, can 
be written as: 


in 

2 K 


( 


1 


A 


r(p,Q)' 


(BA) 
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Appendix C 


C. .1 Two-Dimensional Thermo-elastic Kernels : 


Temperature variation causes thermal stresses in a domain which causes an extra 
term to be introduced into the stress tensor as follows: 


where, the superscript (t) and (e) refer to the thermo-elastic and elastic stresses, 
respectively. Multiplying the stresses by unit outward normal, n/Q), and using 
the traction definition t - a n results in the following expression 


(C2) 

Substituting this thermo-elastic traction and applying the divergence 
theorem into the BIE of Eq. (3.8) results in the following integral equation 

u,(p)+ \T, J (p,Q)u J (Q)dS(Q)= (C.3) 

s 

\u, 1 (p,Q)t J (Q)dS(Q) + ^- f U ^’ q) mdVl3) 

The above equation is valid for both 2-D and 3- applications. For 


Convenience, referring the body force integral as: 

aE cdU.'jfaq) 






dx , 


-mdV(q) 


(C.4) 


Using the expression (3.9) for the two-dimensional kernel Uy the volume 
integral is now a domain integral and can be expressed in terms of distance as: 

= (C5> 

where, A is the area of the enclosed domain. In order to achieve the 
transformation from a domain to a surface integral by using the divergence 
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theorem, some algebraic manipulation is necessary. The domain integral must be 
expressed as a differential operator in V 2 by using the following expression: 

V 2 [r 2 (p,q)\nr(p,q)\ = 4[lnr( M ) + l] (C6) 

therefore, 

V 2 [r 2 ( p,q)lnr(p,q)]-4 


In r(p,q) 


(C-7) 


in: 


Substituting the above integral into domain integral of Eq. (C.5) results 

B ' = ^\ (C ' 8) 
rearranging the differential operator, 

Bi = r(/7 ’ ^)]} 

The above expression is now convenient for use in Green’s second 
identity, which can be stated as: 

[\m~ (C - 10) 

where r(Q) is a line integral. By using the substitution 

A = ~~[r 2 (p,q)\nr(p,q)j (C ‘ ll) 

dx i L • - 1 

and,V 2 ^ = 0, the domain integral of Eq. (C.9) can now be transformed into an 
line integral as follows: 


$ 

+ 

1 ^ 

f 3 

8;r(l -u) p dn 

[dx ; 

a(l+u) ,| 

' 8 j 

8/r(l-u) J| 

dxj ' 


dm) 

dn 


dT(jQ) 
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The above equation contains only surface integrals and that the interior 
point q has now become a surface point Q. For convenience, the thermo-elastic 
surface integrals can be expressed as thermo-elastic kernels as: 

B, = \M,(p,Q)m)dr(Q) + jN,(p,Q)^-dr(Q) (C ' 1 3) 

Therefore, the two-dimensional thermo-elastic kernels can be written as 
follows: 




These kernels can be further expanded to 


M ; 


(l+t>) 


dr{p,Q)dr(p,Q) 




(C.14) 
(C- 15) 

(C.16) 

(C.17) 
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Appendix D 


D.l Implementation Details for Thermal Grooving: 

The curvature value at any node is found out by approximating a circle at 
that node and two of its neighboring points. Thus, the magnitude of the 
curvature, ic is approximately given by: 



where x t c and y t e are the nondimensionalized coordinates of the center of circle, 
given by: 

_ c (y, - 7, .,)[*£■ +3%. -y,)[x?+yf-x U - j£,] 

X ‘ 2 [(* M -^)(7 -y M )-(x, ~x M )(y ,. , -7,)] 

i - (D3) 

y ‘ 2[(? m -x,){y, - 7 )] 

The velocity of groove is obtained by Eq. (4.1) and the Eqs. (4.2) and 
(4.3). (dx/ds) tip can be approximated as {x N -x^)/^ where J N is the non- 

dimensional distance between points “N-l” and “N”. 

The sign of the curvature is found out by drawing a straight line between 
the anterior and the posterior neighbors of the corresponding points as shown in 
Fig. D.l; if this point is located at the right of the line, then a positive sign is 
given to the curvature and vice versa. 

As the factors W, (kTW 4 ) /(D s S s Cly s ) , yJW represent the dimension of 

length, time and stress, respectively, the grooving equation can be normalized on 
the basis of the following relations: 
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Fig. D. 1 : Positive Curvature at node i 

T = (D s S s aY s t)KkTW%s = s/W,x = x!W, (D.4) 

a—a/W, (D.5) 

v„=v„kTW 3 KD s SPr,l (D.6) 

and, 

^00 = crjr/r t , (D.7) 

The non-dimensional form used in the present simulation are given by 

v n =d 2 xlds 2 (D.8) 

cosy/ = y b /(2y s ) (D.9) 

=^i "(2/3)(l ~a)F(dx/ds) lip P.10) 

The non-dimensional normal velocity of atoms at each node, vj, can be 
calculated by 

v, = [(«>+! -(*, +^)/2] (D.U) 

where, I, is the non-dimensional distance between points “i-1” and “i”. 
Consequently, the new coordinates of node points (xj. ’, y f *) are determined by 

(3c ' i \y i ') = (3q +ATv j sin0 n y + ATv i cos0 i ) (D.12) 

9, = tan' 1 (-Ov, -y M ) /(x M -x M )) (D.13) 
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where, (x^y,) are the old coordinates and AT is the non-dimensional time 
increment. The increment of displacement, A Tv i9 is kept below 0.001 at all 
nodes for minimizing the error. The new coordinates at the groove tip are 
obtained by, 

V = Vi '+5V-! ’/t ® V ^ H , 

7 / = 0 "" ' 

An additional node is introduced at the midpoint of the element when the 
distance between the two end nodes exceeds twice the initial spacing. The 
curvature at the additional node is taken to be the average of the curvatures at 
the two neighboring nodes. Repeating the above procedure, the grooving process 
can be fully simulated. This is similar to simulation conducted on the cavity 
growth in the creep by Martinez and Nix [5]. 
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